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SYSTEM FOR COMPUTATIONALLY EFFICIENT ADAPTATION OF 
ACTIVE CONTROL OF SOUND OR VIBRATION 

This application claims priority to U.S. Provisional Application Serial No. 60/271,785, 
Filed February 27, 2001. 

BACKGROUND OF THE INVENTION 

Field of the Invention 

This invention relates generally to improvements in control processes used in active 
control applications, and active control of sound or vibration. More particularly, this invention 
reduces the computations associated with the adaptation process used to tune a controller to 
accommodate system variations by using a more efficient algorithm to implement sound and 
vibration control logic. 

Background Art 

Conventional active control systems consist of a number of sensors that measure the 
ambient variables of interest (e.g. sound or vibration), a number of actuators capable of 
generating an effect on these variables (e.g. by producing sound or vibration), and a computer 
which processes the information received from the sensors and sends commands to the actuators 
so as to reduce the amplitude of the sensor signals. The control algorithm is the scheme by which 
the decisions are made as to what commands to the actuators are appropriate. The amount of 
computations required for the control algorithm is typically proportional to the frequency of the 
noise or vibration. 

Many active noise or vibration control problems, particularly those involving high 
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frequency disturbances, have significant changes in the transfer function between actuator 
commands and sensor response over the system operating regime. Adaptation to these changes is 
required to maintain acceptable performance. The computational requirements associated with 
the adaptation process can be unduly burdensome. Therefore, what is needed is a system that 
reduces computational requirements to implement an adaptation process sufficiently rapidly to 
maintain performance in the presence of a rapidly time-varying system. 

SUMMARY OF THE INVENTION 

The present invention is directed to an apparatus and method for reducing sensed physical 
variables using a more efficient method for updating the transfer function. The method includes 
sensing physical variables and generating control commands at a control rate as a function of the 
sensed physical variables. An estimate of a relationship between the sensed physical variables 
and the control commands is generated, and this estimate is used in generating the control 
commands. At an adaptation rate less than or equal to the control rate, the estimate of the 
relationship is updated based upon a response by the sensed physical variables to the control 
commands. If the control commands are chosen to minimize a quadratic performance metric, 
then the update to the control commands is normalized to maintain constant convergence rates in 
different directions. This normalization factor is inversely dependent on the square of the transfer 
function. To minimize computations, this normalization factor can be updated less often than the 
adaptation rate. This method may be used to reduce vibrations in a vehicle, such as a helicopter. . 

Another embodiment of the present invention is directed to minimizing the computations 
of the control algorithm by updating the quadratic term that the normalization factor depends on, 
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instead of recomputing it when the system estimate is updated. The invention ensures numerical 
stability of this update. . 

Yet another embodiment is directed to directly updating the normalization factor, rather 
than updating the quadratic term on which it depends. The normalization factor can be 
represented as a QR decomposition. The QR factors can be directly updated using a square root 
algorithm. One advantage of this technique is that the normalization factor will always be 
positive definite, that is, that theoretically negative feedback gains are computed as negative 
feedback gains. 

BRIEF DESCRIPTION OF THE FIGURES 

FIG. 1 shows a block diagram of the noise control system of the present invention. 
FIG. 2 shows a vehicle in which the present invention may be used. 

DETAILED DESCRIPTION 

Control systems consist of a number of sensors which measure ambient vibration (or 
sound), actuators capable of generating vibration at the sensor locations, and a computer which 
processes information received from the sensors and sends commands to the actuators which 
generate a vibration field to cancel ambient vibration (generated, for example by a disturbing 
force at the helicopter rotor). The control algorithm is the scheme by which the decisions are 
made as to what the appropriate commands to the actuators are. 

FIG. 1 shows a block diagram 10 of an active control system. The system comprises a 
structure 102, the response of which is to be controlled, sensors 128, filter 1 12, control unit 107 
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and actuators (which could be force generators) 104. A disturbance source 103 produces 
undesired response of the structure 1 02. In a helicopter, for example, the undesired disturbances 
are typically due to vibratory aerodynamic loading of rotor blades, gear clash, or other source of 
vibrational noise. A plurality of sensors 128(a)...(n) (where n is any suitable number) measure 
the ambient variables of interest (e.g. sound or vibration). The sensors (generally 128) are 
typically microphones or accelerometers, or virtually any suitable sensors. Sensors 128 generate 
an electrical signal that corresponds to sensed sound or vibration. The electrical signals are 
transmitted to filter 112 via an associated interconnector 144(a).. .(n) (generally 144). 
Interconnect 1 44 is typically wires or wireless transmission means, as known to those skilled in 
the art. 

Filter 1 1 2 receives the sensed vibration signals from sensors 128 and performs filtering on 
the signals, eliminating information that is not relevant to vibration or sound control. The output 
from the filter 1 12 is transmitted to control unit 107, which includes adaptation circuit 108 and 
controller 106, via interconnector 142. In the present invention, a filter 109 is placed before 
adaptation circuit 108, as will be described below. The controller 106 generates control signals 
that control force generators 1 04(a)... (n). 

A plurality of actuators 1 04(a)... (n) (where n is any suitable number) are used to generate 
a force capable of affecting the sensed variables (e.g. by producing sound or vibration). Force 
generators 1 04(a)... (n) (generally 104) are typically speakers, shakers, or virtually any suitable 
actuators. Actuators 1 04 receive commands from the controller 1 06 via interconnector 1 34 and 
output a force, as shown by lines 1 32(a)... (n) to compensate for the sensed vibration or sound 
produced by vibration or sound source 103. 
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The control unit 1 07 is typically a processing module, such as a microprocessor. Control 
unit 107 stores control algorithms in memory 105, or other suitable memory location. Memory 
105 is, for example, RAM, ROM, DVD, CD, a hard drive, or other electronic, optical, magnetic, 
or any other computer readable medium onto which is stored the control algorithms described 
herein. The control algorithms are the scheme by which the decisions are made as to what 
commands to the actuators 104 are appropriate, including those conceptually performed by the 
controller 107 and adaptation circuit 108. Generally, the mathematical operations described in 
the Background, as modified as described below, are stored in memory 1 05 and performed by the 
control unit 107. One of ordinary skill in the art would be able to suitably program the control 
unit 1 07 to perform the algorithms described herein. The output from the adaptation circuit 108 
can be filtered before being sent to the controller 107. 

For tonal control problems, computations can be performed at an update rate lower than 
the sensor sampling rate as described in co-pending Patent Application entitled "Computationally 
Efficient Means for Active Control of Tonal Sound or Vibration", which is commonly assigned. 
This approach involves demodulating the sensor signals so that the desired information is near 
DC (zero frequency), performing the control computation, and remodulating the control 
commands to obtain the desired output to the actuators. 

The number of sensors is given by n s and the number of force generators is n a . The 
complex harmonic estimator variable that is calculated from the measurements of noise or 
vibration level can be assembled into a vector of length n s denoted z k at each sample time k. The 
control commands generated by the control algorithm can likewise be assembled into a vector of 
length n a denoted u k . The commands sent to the force generators are generated by multiplying 
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the real and imaginary parts of this vector by the cosine and sine of the desired frequency. 

In the narrow bandwidth required for control about each tone, the transfer function 
between force generators and sensors is roughly constant, and thus, the system can be modeled as 
a single quasi-steady complex transfer function matrix, denoted T. This matrix of dimension n s 
by n a describes the relationship between a change in control command and the resulting change in 
the harmonic estimate of the sensor measurements, that is, Az k = T Au k . For notational 
simplicity, define y k = Az k , and v k = Au k . The complex values of the elements of T are 
determined by the physical characteristics of the system (including force generator, or actuator, 
dynamics, the structure and/or acoustic cavity, and anti-aliasing and reconstruction filters) so that 
Tij is the response at the reference frequency of sensor i due to a unit command at the reference 
frequency on actuator j. Many algorithms may be used for making control decisions based on 
this model. For example, one active noise and vibration control (ANVC) approach is described 
below. The control law is derived to minimize a quadratic performance index: 

J = z H W z z + u H W u u + v H W 6u v 
where W z , W u and W 6u are diagonal weighting matrices on the sensor, control inputs, and rate of 
change of control inputs, respectively. A larger control weighting on a force generator will result 
in a control solution with smaller amplitude for that force generator. 

Solving for the control which minimizes J yields: 

u k+1 = u k - Y k (W u u k + T k H W z z k ) 

where 

Yk^CT^WzTk + W. + Wa.)- 1 

Solving for the steady state control (u k+ i = u k ) yields 
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u = -(T H T + W u y 1 T H z 0 

[20] The matrix Y k determines the rate of convergence of different directions in the control 

space, but does not affect the steady state solution. This recursive least-squares (RLS) control 
law attempts to step to the optimum in a single step, and behaves better with a step-size 
multiplier (3<1 . A least means square (LMS) gradient approach would give Y k = I. For poorly 
conditioned T matrices, the equalization of convergence rates for different directions that is 
obtained with the RLS approach is critical. Decreasing the control weighting, W u , increases the 
low frequency gain, and decreasing the weighting on the rate of change of control, W§ U5 increases 
SI the loop cross-over frequency (where frequency refers to the demodulated frequency). 

Sji [21] The performance of this control algorithm is strongly dependent on the accuracy of the 

. estimate of the T matrix. When the values of the T matrix used in the controller do not 

accurately reflect the properties of the controlled system, controller performance can be greatly 
degraded, to the point in some cases of instability. 
[22] An initial estimate for T can be obtained prior to starting the controller by applying 

commands to each actuator and examining the response on each sensor. However, in many 
applications, the T matrix changes during operation. For example, in a helicopter, as the rotor 
rpm varies, the frequency of interest changes, and therefore the T matrix changes. For the gear- 
mesh frequencies, variations of 1 or 2% in the disturbance frequency can result in shifts through 
several structural or acoustic modes, yielding drastic phase and magnitude changes in the T 
matrix, and instability with any fixed-gain controller (z. e., if T M = T k for all k). Other sources of 
variation in T include fuel burn-off, passenger movement, altitude and temperature induced 
changes in the speed of sound, etc. 
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There are several possible methods for performing on-line identification of the T matrix, 
including Kalman filtering, an LMS approach, and normalized LMS. A residual vector can be 
formed as 

E = y - T v 

where no notational distinction is made between the estimated T matrix (available to the control 
algorithm), and the true physical T matrix; all of the control equations are assumed to be 
computed with the best estimate available. The estimated T matrix is updated according to: 
Tk+i = Tk + EK H 

The different estimation schemes differ in how the gain matrix K is selected. The 
Kalman filter gain K is based on the covariance of the error between the true T matrix and the 
estimated T matrix. This covariance is given by the matrix P where P 0 = I, and 

M = P k +Q 

K = Mv/(R + v H Mv) 

Pk+i^M-K^M 

and the matrix Q is a diagonal matrix with the same dimension as the number of force generators, 
and typically with all diagonal elements equal. The scalar R can be set equal to one with no loss 
in generality provided that both Q and R are constant in time. The normalized LMS approach is 
simpler, with the gain matrix K given by: 
K = Qv/(l+v H Qv) 

The computational burden associated with updating T k is roughly 2n a n s (using the 
normalized LMS gain rather than the Kalman filter gain). This is not overly burdensome, and 
cannot be readily avoided. However, the update equation for u k +i requires not only T k , but the 
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triple product T k H W z T k and the inverse (T k H W z T k + W u + W 5 u) _1 . These two steps are 
computationally intensive, but potentially amenable to some straightforward investigation. First, 
the inverse need not be computed directly. Since Y k 4 = (T k H W z T k + W u + W 5u ) is Hermitian, the 
required product can be obtained by first computing the Cholesky decomposition, from which the 
required product can be obtained by back-substitution. The Cholesky decomposition requires 
roughly n a 3 /6 floating point operations (flops), plus computations on the order of n a 2 . Another 
significant modification that appears to be straightforward is to propagate X k = T k H W z T k , rather 
than computing the matrix multiplication at each step . Given that T has a rank one update, T k+! = 
T k + EK H , then X k+i satisfies 

X k+ i = X k + (T k H W z E)K H +K(T k H W z E) H +K(E H W z E)K H 
[-26] However, without further modification, this equation is numerically unstable and cannot 

m . ■ 

be implemented. Random numerical errors due to round-off or truncation that are introduced at 
|l each step accumulate until eventually, X k diverges from T k H W z T k , potentially leading to 

fui 

\. i instability of the overall algorithm. 

f : 'l 

|y [27] Without modifications, the computations of the overall algorithm remain significant, and 

for many applications, the resulting burden is unacceptable. An algorithm is desired that gives 
equivalent performance, with much lower computation. 
[28] One embodiment of the present invention is directed to reducing the computational 

burden. The primary difficulty with the baseline algorithm for noise control is the computational 
burden. This is driven by the computation of T H T, and by the solution of the equation for u. 
Assume that W u , W 6u , W z and Q are all diagonal matrices. If the matrix-multiplication is 
computed directly, and a Cholesky decomposition used to solve for u, then the computational 
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burden of the algorithm in flops is roughly n a 3 /6 + n a 2 n s + 3n a 2 + 3n a n s , ignoring vector 
computations which are linear in n a or n s . As noted in the algorithm derivation, the matrix Y k 
affects only the convergence rate, and not the converged solution. Therefore, it does not need to 
be updated at the same rate as the control and adaptation. Splitting the computation of the 
Cholesky decomposition over several control iterations reduces the computations per iteration. 
For example, the Cholesky decomposition can be split over 4 steps. Performing all of the 
adaptation at a lower rate is also possible. However, noting that the two different uses of the 
estimated T matrix (Le. for computing the gradient, and for normalizing the directions) result in 
different demands on the accuracy of T leads to better use of the available computation. The 
matrices W u and W Su can be time varying, but can only be changed during an iteration when the 
Cholesky decomposition is updated (that is, the W u used to compute a must be the same as that 
used to compute the Cholesky factors). 

Another embodiment of the invention is directed to using the update equation for X. 
Since numerical errors will always be introduced at every step, over time, X k will gradually 
diverge from T k H T k . (The dynamics associated with the propagation of numerical error in the 
above equation are neutrally stable.) To prevent this divergence, the update equation for X can 
be modified so that it depends on X itself. The form of the above update equation will guarantee 
that X is positive definite and Hermitian, and any modification must maintain this behavior. 
Noting that T H W Z E = T H W z y - T H W z Tv, then define instead E x = T H W z y - X v and substitute this 
into the previous update equation for X. The resulting equation will still guarantee that X k+ i is 
positive definite and Hermitian, and X will still satisfy X = T H W Z T except for numerical errors. 
However, an analysis of the error propagation reveals that the error behavior is now strictly 
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stable, and thus cannot accumulate indefinitely. 

Another embodiment of the present invention is a more efficient computation for a 
control update algorithm. The definition of E x above involves T k H W z y = T k H W z z k -T k H W z z k _t. 
Since the control update equation already computes F k = T k H W z z k , then E x can be computed as: 

E x = F k -F k _ r F c -Xv 

where the correction term F c is given by ¥ c =K k .{E k .i U W z z k . h This computation involves only 
vector computations, and is thus efficient. 

2 * 

The update equation for X k +i involves 3 terms, each corresponding to n 12 computations 
accounting for symmetry. However, these terms can be grouped to form 2 rank 1 updates, rather 
than 3. Modifying the definition of E x gives us: 

E x = F k -F k4 -F c -Xv+(E H W z E)K/2 

X k+ i=X k +E x K H +KEx H 

The above equations assume that W z is diagonal and constant. However, if W z is allowed 
to be time-varying, then the update equations for X must change. If complete freedom is allowed 
in the time variation, then no computational simplifications from the above steps can be applied. 
However, if one permits only a single element of W 2 to change at each iteration, then the change 
in X can be computed via a computationally efficient rank one update. If the kth element of W z 
increases by (AW z ) k , then the modification to X can be computed as follows, where T k refers to 
the k th row of the T matrix: 

X new = Xoi d + (AW 2 ) k T k H T k 

Examining the behavior of the adaptation, the diagonal elements of the covariance are 
most important, and the off-diagonal elements have little impact on performance. Making the 
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covariance a real vector consisting of only the diagonals saves 2n a 2 operations. Further 
simplifications to eliminate the time-varying covariance P results in an equation identical to the 
normalized LMS approach described previously. 

Incorporating all of the above modifications results in an algorithm with roughly 7n 
operations per step; an improvement of roughly a factor of 6 over the original algorithm, with 
almost no change in the behavior of the algorithm. To summarize, the new equations are as 
follows: 

F k = T k H W z z k ; 

S H S = Chol(X k + W U + W 5u ) (every 4 iterations); 

Uk+1 = Uk . (S H S)" 1 (W u u k + F k ) (the product is computed via back-substitution); 

v=Au; 

y=Az; 

E=y-T k v; 

K = Q/(l+v H Qv); 

Ticfi = T k +EK H ; 

F c = K k .,E k .i H W z z k . 1 ; 

E x = F k - F k ., - F c - Xv + (E H W z E)K/2; 

X k+1 = X k + ExK H + KE X H ; and 

X ne w = Xoi d +(AW z ) k T k H T k . 

Ignoring vector and scalar operations, the total computational burden associated with the 
current algorithm is: 

Control update: 1 matrix-vector multiply (n a n s ) 
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Cholesky back-substitution (na 2 ) 



Cholesky decomposition: 



n a 3 /6, split over 4 steps 



Residual calculation: 



1 matrix- vector multiply (n a n s ) 



Adaptation filter gain: 



vector operations only 



Update of T estimate: 



1 vector outer product (n a n s ) 



Computation of Ex: 



1 matrix- vector multiply (n a 2 ) 



Computation of X: 



2 symmetric outer products (n a 2 ) 



sym. outer product for variable W z (n a 2 /2) 

Another embodiment of the present invention is directed to improving the efficiency of 
calculations by using a square-root algorithm that enables a controller 106 to achieve the same 
attenuation of a physical variable, such as noise, sound or vibration while using less expensive 
computer hardware. Alternatively, the same computer hardware can be used to control 
approximately twice as many modes of vibration or sound. This algorithm achieves the same net 
computation precision as algorithms for quasi-steady control logic, but with computer hardware 
that is only half as precise in each operation. For example, if double precision, floating point 
arithmetic is required for a particular control algorithm, this algorithm would only require single 
precision arithmetic. Since single precision operations are much faster, the same controller can 
be implemented with a slower, less expensive computer. The algorithm described herein allows 
lower cost active noise and vibration control systems. 

In addition to doubling the precision, the algorithm described herein is an inherently more 
stable implementation. In conventional algorithms, numerical errors can cause modes that are 
theoretically stable to become unstable. For these modes, the numerical errors cause slightly 
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stable negative feedback gains to be computed as slightly positive feedback gains and, thus, they 
become unstable. Due to the nature of the numerical method in the square root algorithm, 
theoretically negative feedback gains are computed as negative feedback gains despite numerical 
errors. 

Most active controllers of sound and/or vibration are based on quasi-steady control logic. 
That is the source of the sound and vibration is a persistent excitation of one or more discrete 
frequencies that vary relatively slowly. The amplitudes and phases of the discrete frequencies 
take one or more seconds to change significantly. The algorithm described herein applies to 
quasi-steady control logic. 

Quasi-Steady Control Logic 
Quasi-steady control logic refers to optimal control logic for multi-variable systems assumed to 
have transfer functions that do not vary within the frequency range that needs to be controlled. 
Quasi-steady control logic is commonly used in sound and vibration control because the transfer 
functions relating actuator inputs to microphone or accelerometer outputs do not vary 
significantly in a narrow frequency band about the frequency of one of the discrete frequency 
disturbances. If there are multiple discrete tones that need to be attenuated, the controller would 
use a separate control logic for each. For each tone, the system is modeled by a transfer function 
that consists of a matrix of constant gains. For convenience, the m inputs, u n , and p outputs, y n , 
are modeled with separate real and imaginary parts and thus the p x m transfer function matrix, T, 
consists of real numbers. Alternatively, complex gains could be used. 

The optimal control problem is to minimize the performance index, J , at time n through 
selection of a perturbation, Au n to the control signal, where: 
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Jn = 0.5*(Yn T *yn + Au n T *W*Au n ); 
y n = TAu n + y n .i;and 
u n = u n .i + Au n . 

W is a real and positive semi-definite m x m control effort weighting matrix. The optimal 
control is that which causes: 

5Jn/5Au n = (Sy n /SAu n ) T *y n + W*Au n = 0. 
This implies the optimal control is: 
Aun^-CT^T + Wy^T 1 ^!. 

In noise and vibration control the control logic is made adaptive by estimating the values 
of T. As discussed herein, T refers to the estimate of the transfer function matrix. Assuming that 
each element of the transfer function is a Brownian random variable, the minimum variance 
estimate of it at time n+1, T n +i, is: 

T n +i = T n + E n *L T , 

where E n = y n - yp n are the innovations, yp n = T n *Au n + y n -i, is the predicted value of y at time n, 
and L is a m x 1 matrix of constant gains. This type of estimator is a Kalman filter. 
In summary, the adaptive quasi-steady control logic is: 

T n = T n .i+(yn-yPn)*L T , 

Au n = -(T n *T n +W)- 1 *T T n *y n (1) 
yPn+i=yn + T T n *Au n 

U n =U n + Au n 

Formulation as a QR Problem 

The control logic can be reformulated in terms of a matrix decomposition into the product 
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of a orthonormal matrix, Q, and a upper triangular matrix, R. This is called a QR decomposition. 
The symmetric, positive definite m x m matrix, Y n will be decomposed and propagated via a 
square root method where: 

YnKW+T/Tn)" 1 

Propagating Y n 

Y n can be propagated using the following recursive relationship. Combining the 
definition of Y n and the Kalman filter update for T n yields: 

Y n 4 = Y n .f l + LE n T T n .i+T n .i T E n L T +LE n T E n L T ? 
which can be more compactly expressed as: 

Y n _1 = Y n .f 1 + c n p V-b n -iP~Vi T , 
using the definitions: 

C n = T n EnJ 

d n .i=T n .i T E n ;and 
p 2 = E n T E n . 

Collecting the time n terms of the Y propagation equation into the left hand side, 
inverting both sides of the resulting equation, and using the matrix inversion lemma yields the Y 
propagation equation: 

Y n +Y n Cnr n 2 c n T Y n = Y n . 1 +Y n -id n .iv n _ 1 2 d n .i T Y n . b (2) 
where r n and v n -i are defined as: 

tn 2 = (p 2 -c n T Y n c n )" 1 ;and 

Vn.i'^Cp'-dn-l^ldn-l)' 1 . 

To present this as a QR decomposition, each term must be expressed in the quadratic 
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form c T c, where c is real. Since Y n and Y n +i are positive semi-definite and symmetric, real upper 
triangular matrices, R n and R n .i can be defined such that: 
R n T Rn = Y„;and 

Rn-l T Rn-l = Y n .i 

These are known as a Cholesky decompositions. Putting the remaining terms in quadratic 
form only requires that, r n and v„_i, be real. Using the definitions of Y, c, and r, 
r n 2 = p 2 -E n T T n (W+T n T T n y VE n 
= En T 0[-T n (W +T n T T n )" 1 T n T )E n 

= En^-TnW^I+T^TnW-V^^En 

=E/((I+T n W- 1 T T n )- I (I+T n W- 1 T\)-T n W- 1 T T n a+T n W- 1 T T n )- 1 )E n 
=E„ T (I+T n W 1 T T n )- 1 E n . 

This result is positive because the matrix within the parenthesis is symmetric and positive 
definite. Thus r n will be real. v„.i can be shown to be real following the same procedure. 
The Y propagation equation can be put in the following quadratic form: 









r - 
* 


L *» . 









I *-> J" 



This can be put in the form of QR decomposition by adding an appropriate column vector 
as follows: 



[o R n J y It o vi J U 'J. 



(3) 



where Q is an orthonormal matrix. If each side of Equation (3) is multiplied on its left by its 
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transpose, the equation above is one if the results. However, Equation (3) allows the following 
algorithm to be used for the propagation. Starting with the upper triangular matrix on the right 
hand side of Equation (3) from time n-1 it is converted to the first matrix on the left hand side of 
the time n equation replacing the first row with the terms shown. This is how the new 
information inherent in the measurement y n is entered into the square root propagation. Next, it 
is multiplied on the right by the matrix containing L. 
[52] Finally, a series of orthonormal row operations are performed on the resultant matrix to 

produce an upper triangular matrix. These row operations can be collected into the form of an 
IS orthonormal matrix, Q T , pre-multiplying the matrix. This final operation is termed a QR 

decomposition. The resulting upper triangular matrix has the form of the time n-1 result, but 
with its values updated to time n. Q does not need to be actually formed. Instead of 
propagating Y, its square root, R, is propagated instead. For this reason the numerical precision 
needed to propagate Y in a computer implementation is reduced by approximately half. 
\| [53] The control logic contains the term Y n T„ T yn. This can be put in terms of one of the results of the 
QR decomposition, saving some computations. 

Y n Tjy n = Y n T^(E n +yp n ) = Y n T^E n + Y n (T^ + LE^)yp n = Yrfn-YnCWAun-i-UBSyPn) 

using 

TLypn = (I-TLTn.iCW+T^T.^r^Tliyn-i = WCW+TLTnO'TLyn^ -WAu^ 
[54] The remaining control algorithm, including the Kalman filter is: 

Au n = -Y n r n +Y n (WA n _i-LEn T yp n ) 

Tn-Tn.j+En*!/, (4) 

yPn+i=yn+T n Au n . 
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Equations (3) and (4) are the control logic of Equation (1) reformulated as a QR 
decomposition. 

These QR equations can be confirmed by multiplying each side of the equation on the left 
with their respective transpose matrix. This yields a block symmetric matrix equation with the Y 
propagation equation, Equation 2, appearing in the lower right block. It remains to show that the 
off-diagonal and upper left blocks are consistent with Equation 2. 

The off-diagonal submatrix from the right hand side is 

(1 +L T Y n _id n . 1 )viidliY n .i+L T Y n .i 

= viidl 1 Y n .i+p" 2 (cI-dLi)(Y n . 1 +Y n . 1 d n .iv^ 1 dIiY n .i) 
where E n was expressed in terms of c n and d n _i. Factoring the above into c n and d n _i, components 
yields 

p- 2 cI(Y n 4+Y n 4d^ 

The second term is zero. Substituting in Equation (2) into the first term yields 

p- 2 cI*(Y n +Y n c n ^Y n ) = p- 2 (l+c^Y n *c n *r 2 )*c^Y n = r 2 *d*Y n . 

Which is the off-diagonal term on the left hand side of Equation (3). 

The upper left submatrix from the right hand side of the QR formulation is 

(l+L T Y n .ib n -i)qn-i(l+bIiY n _iL)+LViL 

Substituting in the relation to d n .i, and c n for the post multiplication by L , and factoring 

yields 

((l+L T Y n _ 1 d n4 )vi 1 di 1 Y n _ 1 +L T Y n _ 1 )c n p- 2 
+(l+L T Y n . 1 d n . 1 )vii(p 2 -dliY n . 1 d n . 1 )p- 2 -L T Y n . 1 d n .ip- 2 

The term in the outside parentheses is the off-diagonal term. Substituting in the off- 
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diagonal result and using the definition of q£ twice results in 
(l+^Y n c n )p- 2 = r^ 

Which is the upper left submatrix on the left hand side of Equation (3). 
Modified Givens Method 

Any matrix can be decomposed into an orthonormal, matrix, Q, pre-multiplying an upper 
triangular matrix, R. In Equation (3) the (m+1) x (m+1) matrix to be decomposed: 



L 0 R *-i 





"1 0' 


* 


L. I 



is almost in upper triangular form. The only exception is that the first column has nonzero 
entries. A matrix in this form can be decomposed into Q and R with far fewer computations than 
required for a general matrix. The following approach modifies the known Given's method of 
QR decomposition for a general matrix to exploit the sparsity of the lower triangular portion of 
the above matrix. Decomposition is accomplished by choosing Q to consist of a sequence of 
Given's transformations. Each Given's transform produces one zero in the matrix, by operating 
on two matrix rows with a Given's rotation. Each Given's transform has the form 
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0 



where G, ~ 



-b a 



Then 
G 



[x — 



x a x •• 
x b x 



x x Va 2 +6 2 x - x 
x .« x 0 x x 



The sequence of Given's rotations consists of a reverse pass sequence followed by 
forward pass sequence. The first Given's rotation operates on the last two rows of the matrix to 
make the last row of the first column zero. The next in the reverse sequence operates on the m-1 
and m rows to make the m row of the first column zero, and so on until the 3rd row of the first 
column is zero. The result of this backward sequence of orthonormal transformations is a matrix 
with zeros in the first column as needed, but with nonzero entries along the sub-diagonal below 
the main diagonal. The forward sequence puts these sub-diagonal terms back to zero without 
disturbing the zeros in the first column. 

The first Given's rotation of the forward sequence operates on the first two rows of the 
matrix to make the second row of the first column zero, the next operates on the 2nd and 3rd 
rows to make the 3rd row of the 2nd column zero, and so on until the last row of the second last 
column is zero. The resulting matrix is now in upper triangular form and therefore it is 

\ rc T n Y' 

it u n ft 

0 R„ 
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Note that the orthonormal matrix, Q, does not need to be explicitly computed. The 
number of computer operations required varies with the number of sensors, p, and the number of 
actuators, m. In estimating the number of computer operations only square root operations and 
multiplications and divisions, termed an op, will be counted. Multiplications by zero do not have 
to be done and are not counted. It takes four multiplication's and one square root to determine 
each Givens transformation. Performing the reverse sequence transformation on the j and j+1 
rows requires 2 + 4*(m-j+l) ops, for a total of 10 + 4*(m-j) plus one sqrt. In the reverse 
sequence, this set of operations needs to repeated for j = m, m-1, 2. Thus, the reverse 
sequence requires 2m 2 +4m-6 ops and m-1 square roots. Similarly, the forward sequence requires 
2m 2 +6m ops and m square roots. Thus, the Given's method of QR decomposition for the spare 
matrix requires 4m 2 +10m-6 ops and 2m- 1 sqrts. 

Numerical Stability 

Theoretically, the matrix Y has all positive singular values. However, numerical errors in 
directly computing can result in small positive singular values becoming small negative singular 
values. This might make a theoretically stable sound and vibration control system unstable. The 
square-root method avoids this potential problem by not forming Y but using its square root 
instead. In spite of numerical errors the square root matrix, R, will only contain real values. 
Thus, R T R can have only positive singular values. 

Algorithm and Operation Count 

The algorithm for the n th time point is: 



Operation Sequence 



Op Count 
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E= (yn-ypn) 


0 


p 2 =E T E 


p 


d=T T E 


mp 


R*d 


(m 2 +m)/2 (since R is upper 
triangular) 


v=sqrt(p 2 ~ (R*d) T (R*d) ) 


m+1 sqrt 


R T (R*d) *v 


(m 2 +m) /2+m 


v n +{Y*d*v)n T *L 


m 


Rn*L 


( (m 2 +m)/2 


m+1 order QR 


4m 2 +10m-6 ops and 2m- 1 sqrt 


U n =U n -l- Y n r n +R T R (W_U n -l- 

LE n T YPn) 


m 2 +4m+p 


Tn=T n+ i+E n *L T 


m*p 


ypn + l=Yn+Tn_Un 


m*p 


total 


2m sqrts plus 

(6.5 m 2 +3 m*p+17.5 m + 2p-6) ops 



input data: y„ 
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in memory from n-1 calculations: S n , yPn, 

(q*Z*b) n ,u n . 

constants: L, r, W" 1 (W is assumed to be a diagonal matrix). 

The square root method requires fewer computer operations than other algorithms 
implementing the adaptive quasi-steady vibration and/or noise control logic. The logic, described 
in Equation (1), is repeated here for convenience. 

Au^TVTn+WX^TVyn 

yPn+i=yn+T n *Au n 

u n =u n +Au n 

Simply executing this control logic as shown requires 3m*p+m 2 operations in addition to 
the operations required for forming the matrix inverse. Other than the square root method 
disclosed here, there is no known method for forming the required inverse that uses as few as 5.5 
m 2 + 17.5m + 2p- 6 ops. 

Alternate Formulation 

By substituting TW" 1/2 for T T , W' 1/2 L for E n , E n for L, Z for Y and S for R an alternate 
form of QR formulation can be determined. In the alternate propagates the pxp matrix: 

Z n =(I+T n W 1 T n T )" 1 . 
Using Zn 

Z„ can be used to compute both Au n and yp n . The derivation of the corresponding 
relations, will use the matrix equalities: 
Y(M-XY)" 1 = (I+YX)-'Y, 

anda+vy^i-a+vyV 
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which can be verified by multiplying through by the respective inverted matrices. Using these 
equalities 

Zn^I+TnW-V) 4 = [I-TnWV (I+T n W" V) 4 H-Tn (T^n+WyV 

Comparing this to the control logic above shows that 
yPn-H = Z n y n 

The control Au n+ i can also be expressed in terms of Znt 

AUn+i-W-VZnyn. 

This can be verified using the above matrix equalities once again. 

-WV Z n y n = -W'V (I+TnW^n^Vn = "(TVlTn+l+^TnViyn^Un+i 

Applying the substitutions listed above to the Y propagation equations yields the Z 
propagation equations 

2 T 2 T 

Zn+Znbnqn b n Z n = Zn-l+Zn-lb n .iq n -l b n -i Zn-i, 

using the definitions 

q„ 2 = (r^VZnbn)- 1 , bn^W'L, and r 2 = L J W l L 
Then the dual QR formulation is 



q a q x b T H Z„ 
0 S„ 



i 0 



where Z„-i = SViS n -i, yp n = Z n _]y n .i, and E n = y„- yp n , 

ypn+i = S n S n y n 
T n = T n _i + E n *L T , 
Au n =-W 1 *T n T *yp n+1 . 

The alternative form has the advantage that after the substitutions v n =r n , d n =c m and r is 
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constant. Therefore the computations shown in the table rows 2 through 6 do not need to be 
performed. It has the disadvantage that the QR decomposition is on a p+1 square matrix rather 
than the normally smaller m+1. The op count for the alternative formulation is found by 
switching the roles of m and p in the remainder of the table: 5.5p 2 + 2mp + 12.5p +m -6 ops. 
Generally, this form only has an advantage in operation count if p < 1 .1 8*m. 
[80] Adaptive quasi-steady vibration and/or noise control with square-root filtering is 

extremely attractive for implementation. The square root algorithm can provide a desired level of 
computation performance with significantly less computer power. It is also more numerically 
stable. 

*3 [81 ] FIG. 2 shows a perspective view 20 of a vehicle 1 1 8 in which the present invention can be 

fi used. Vehicle 118, which is typically a helicopter, has rotor blades 119 (a)...(d). Gearbox 

HI 

housing 1 10 is mounted at an upper portion of vehicle 118. Gearbox mounting feet 140 (a)...(c) 
JFl (generally 140) provide a mechanism for affixing gearbox housing 1 10 to vehicle airframe 142. 

Sj Sensors 128(a) through (d) (generally 128) are used to sense acoustic vibration produced by the 

fyt vehicle, which can be from the rotorblades 1 1 9 or the gearbox housing 110. Although only four 

sensors are shown, there are typically any suitable number of sensors necessary to provide 
sufficient feedback to the controller (not shown). The sensors 128 may be mounted in the vehicle 
cabin, on the gearbox mounting feet 140, or to the airframe 142, or to another location on the 
vehicle 118 that enables vehicle vibrations or acoustic noise to be sensed. Sensors 128 are 
typically microphones, accelerometers or other sensing devices that are capable of sensing 
vibration produced by gear clash from the gearbox 1 10 and generating a signal as a function of 
the sensed vibration. These sensors generate electrical signals (voltages) that are proportional to 
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the local noise or vibration. 

In accordance with the provisions of the patent statutes and jurisprudence, exemplary 
configurations described above are considered to represent a preferred embodiment of the 
invention. However, it should be noted that the invention can be practiced otherwise than as 
specifically illustrated and described without departing from its spirit or scope. Alphanumeric 
identifiers for steps in the method claims are for ease of reference by dependent claims, and do 
not indicate a required sequence, unless otherwise indicated. 
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